/*
 * TopolPDB.cpp
 *
 *  Created on: Aug 24, 2013
 *      Author: marchi
 */

#include "TopolPDB.h"

namespace Topol_NS {
struct op_sortPDB{
	bool operator()(string x, string y){
		stringstream ssx(x.substr(21,5));
		stringstream ssy(y.substr(21,5));
		int nx, ny;
		ssx>>nx; ssy>>ny;
		return (nx < ny);};
};

vector<string> split(const string & s){
	stringstream iss(s);
	vector<string> tokens;
	copy(std::istream_iterator<string>(iss),
			std::istream_iterator<string>(),
			std::back_inserter<vector<string> >(tokens));
	for(unsigned int o=0;o<tokens.size();o++)
		tokens[o].erase(remove_if(tokens[o].begin(),tokens[o].end(),::isspace),tokens[o].end());
	return tokens;
}

void TopolPDB::operator()(const vector<string> & data0 ){
	map<string,vector<string> > mapps;
	vector<string> data=vector<string>(data0.size(),string());
	int ia=0;
	for(size_t i=0;i<data0.size();i++){
		if(data0[i].substr(0,7).find("ATOM") == 0 || data0[i].substr(0,7).find("HETATM") == 0 ) {
			data[ia]=data0[i];
			ia++;
		}
	}
	data.resize(ia);
	op_ResNo tmp;

	std::transform(data.begin(),data.end(),data.begin(),tmp);

	int nresMax;
	std::stringstream(data[data.size()-1].substr(22,5))>> nresMax; // Residue number
	vector<vector<string> > vstring=vector<vector<string> > (nresMax);
	for(size_t i=0;i<data.size();i++){
		int nres;
		std::stringstream(data[i].substr(22,5))>>nres; // Residue number
		vstring[nres-1].push_back(data[i]);
	}
	size_t ntot=0;
	for(size_t n=0; n< vstring.size(); n++){
		ntot+=vstring[n].size();
	}
	for(size_t n=0; n< vstring.size(); n++){
		vector<PDBdata> tmp=SplitResidue(vstring[n]);
		PDB.insert(PDB.end(),tmp.begin(),tmp.end());
	}
}
vector<PDBdata> TopolPDB::SplitResidue(vector<string> & v0){
	string sub2=v0[0].substr(17,4); // Residue name
	sub2.erase(remove_if(sub2.begin(),sub2.end(),::isspace),sub2.end());

	Mymapmap::iterator it=ResSplit.begin();
	if(!ResSplit.empty()) {
		do{
			if(sub2.find(it->first) != string::npos) break;
		} while(++it != ResSplit.end());
	}

	vector<int> atmindx=vector<int>(v0.size());
	vector<int> resindx=vector<int>(v0.size());
	vector<string> resname=vector<string>(v0.size());
	if(it != ResSplit.end()) {
		map<string,int> original;
		for(size_t o=0; o < v0.size();o++){
			string sub1=v0[o].substr(11,5); // Atom name
			sub1.erase(remove_if(sub1.begin(),sub1.end(),::isspace),sub1.end());
			original[sub1]=o;
		}
		Mymap & ResSplit0=it->second;
		Mymap::iterator ip=ResSplit0.begin();
		int na=0;
		int pp=0;
		for(;ip != ResSplit0.end(); ip++){
			for(size_t o=0;o < ip->second.size(); o++){
				atmindx[pp]=original[ip->second[o]];
				resindx[atmindx[pp]]=na+offsetRes;
				resname[atmindx[pp]]=ip->first;
				pp++;
			}
			na++;
		}
		offsetRes+=na;
	} else{
		for(size_t o=0; o < v0.size();o++){
			atmindx[o]=o;
			resindx[o]=offsetRes;
			resname[o]=sub2;
		}
		offsetRes++;
	}

	vector<PDBdata> y=vector<PDBdata>(v0.size());
	vector<PDBdata> yout=vector<PDBdata>(v0.size());
	for(size_t o=0;o < v0.size(); o++){
		string sub1=v0[o].substr(11,5); // Atom name
		sub1.erase(remove_if(sub1.begin(),sub1.end(),::isspace),sub1.end());
		y[o].atn=sub1;
		std::stringstream ss(v0[o].substr(31,24)); // coordinates
		ss>> y[o].x;
		// And now let us do the subunit stuff
		float Mya;
		std::stringstream(v0[o].substr(61,6))>>y[o].domd;
		y[o].domn=v0[o].substr(21,1);
	}
	for(size_t o=0;o < y.size(); o++){
		int oo=atmindx[o];
		yout[o].at=atmindx[o]+offset;
		yout[o].res=resindx[o];
		yout[o].resn=resname[o];
		yout[o].atn=y[o].atn;
		yout[o].x=y[o].x;
		yout[o].domd=y[o].domd;
		yout[o].domn=y[o].domn;
		yout[o].first=v0[o].substr(0,61);
		yout[o].last=v0[o].substr(66,14);
	}
	offset+=v0.size();
	return yout;
}

TopolPDB::~TopolPDB() {
	// TODO Auto-generated destructor stub
}
std::ostream & operator<<(std::ostream & fout, TopolPDB & y){
	vector<PDBdata> & PDB=y.PDB;

	vector<string> tmp=vector<string>(PDB.size());

	for(size_t o=0; o < PDB.size();o++){
		stringstream ss;
		ss <<std::setw(4) << std::left << PDB[o].resn;
		tmp[o]=PDB[o].first.replace(17,5,ss.str());
		ss.str(string());
		ss << std::setw(5) << std::right << PDB[o].res;
		tmp[o].replace(21,5,ss.str());
		tmp[o].append(" 1.00");
		tmp[o].append(PDB[o].last);
	}
	std::sort(tmp.begin(),tmp.end(),op_sortPDB());
	for(size_t o=0; o < tmp.size();o++){
		fout << tmp[o] <<endl;
	}
	return fout;
}

std::ifstream & operator>>(std::ifstream & fin, TopolPDB & y){
	vector<string> data;
	for(string str;getline(fin,str);){
		data.push_back(str);
	}
	vector<string>::iterator it=data.begin();
	do {
		vector<string> tokens;
		tokens=split(*it);

		string Tag0,Tag;
		if(tokens.size() ==1){
			Mymap ResSplit0;
			if(tokens[0].substr(0,1) == "$") Tag0=tokens[0].substr(1,tokens[0].size()-1);
			it++;
			int m=0;
			stringstream ss;
			ss<<m++;
			do {
				Tag=Tag0+ss.str();
				tokens=split(*it);
				for(size_t n=0;n < tokens.size();n++){
					if(tokens[n].find("next") == string::npos)
						ResSplit0[Tag].push_back(tokens[n]);
					else {ss.str(string());ss<< m++; Tag=Tag0+ss.str();}
				}

			} while((*(++it)).find("end") == string::npos );
			y.ResSplit[Tag0]=ResSplit0;
		}

	} while(++it != data.end());
	return fin;
}

} /* namespace Topol_NS */
